function [alp,bet,mode]=msdbetab(mu,sd,a,b,flag_plot) 
% ================================================================
% function [alp,bet,mode]=msdbetab(mu,sd,a,b)
%
% Given y in [a,b], assume it has a density with mean MU and std SD 
% Model as a Beta(alph,bet) 
% The transformation is 
% y=a+(b-a)*x  when x~Beta(alph,bet) 
% If MU is not given then choose the midpoint 
% 
% Output 
% ALP       First location parameter of Beta
% BET       2nd   location parameter of Beta
% MODE      
% Alejandro Justiniano August 7 2008 
% ================================================================
if nargin < 5;flag_plot=0;end; 
if b <= a; error('A must be smaller than b');end 
if isempty(mu) 
    disp('Choosing MU=midpoint') 
    mu=0.5*(a+b); 
end 
if mu < a; error('MU must be > a');end 
if mu > b; error('MU must be < b');end 
if sd < eps;error('SD must be positive');end 
% FIND ALPHA and BETA 
munew=(mu-a)/(b-a); 
sdnew=sd/((b-a)); 
bet=munew*( (1-munew)^2 )/(sdnew^2)-(1-munew );
alp=munew*bet/(1-munew);
% Density is not well behaved if ALP <= 1 || BET <=1 
% In this case, mode is usually A,B or Bimodal 
if bet <=1; disp('Warning: BETA is <=1');end 
if alp <=1; disp('Warning: ALP  is <=1');end 
if flag_plot==1
    x=linspace(a+eps,b-eps,10000)';
    dens=(betapdf(((x-a)/(b-a)),alp,bet))/(b-a);
    mode=a+(alp-1)*(b-a)/(alp+bet-2);
    plot(x,dens);
    vline(mode);
    xlim([x(1) x(end)]);
    str=['Mean=',num2str(mu),'   Std=',num2str(sd),'   Alpha=',num2str(alp),'   Beta=',num2str(bet),'    Mode=',num2str(mode)];
    title(str,'FontSize',11);
end
% Plot New Density 
% cons=gamma(alp+bet)/(gamma(alp)*gamma(bet)*(b-a)^(alp+bet-1));
% d1=(x-a).^(alp-1); 
% d2=(b-x).^(bet-1); 
% dens1=cons*(d1.*d2); 
% plot( dens2./dens1 ); 